#| Given linear eqns: 
                      x = G x + b, with det(E-G) != 0,
we want to get the unique solution of the eqns by iteration:
                      x_{n+1} = G x_n + b

We know if the radius of spectrum rho(G) > 1, it will not converge.
In this note we consider that after the Aitken transformations, how about the iteration.

|#

(import srfi-41)
(import list-comprehensions)
(load "matrix-vector-operations.scm")
(load "iterative-methods-use-Stream.scm")

;; nextstate
(define (nextstate G b x)
  (vector-map +
              (matrix-*-vector G x)
              b))
(define (multi-nextstate G b x n)
  (let iter ([count 0] [result x])
    (if (= count n)
        result
        (iter (+ count 1) (nextstate G b result)))))
;; stream
(define (general-stream-from G b x)
  (stream-cons x
               (general-stream-from G b (nextstate G b x))))
(define (multi-general-stream-from G b x n)
  (stream-cons x
               (multi-general-stream-from G b (multi-nextstate G b x n) n)))
;; Aitken acceleration
(define (Aitken-general-strm G b x0)
  (Aitken-acceleration (general-stream-from G b x0)))
(define (multi-Aitken-general-strm G b x0 n)
  (Aitken-acceleration (multi-general-stream-from G b x0 n)))

;; super acceleration
(define (super-general-strm G b x0)
  (super-accleration Aitken-acceleration (general-stream-from G b x0)))
(define (super-multi-general-strm G b x0 n)
  (super-accleration Aitken-acceleration (multi-general-stream-from G b x0 n)))

;;-------- test
(define (Jordan c n)
  (vector-map (lambda (i)
                (vector-map (lambda (j)
                              (cond [(= j i) c]
                                    [(= j (+ i 1)) 1]
                                    [else 0]))
                            (list->vector (range n))))
              (list->vector (range n))))
#| Note:
(i) If det(E-G) != 0, where G is Jordan, it must c != 1;
(ii) If |c|<1, the iteration converges, while |c|>1, it fails;
(iii) The calculation below shows that, Aitken-acceleration can not save the fath of iteration.
;; in some cases, it converges if c<-1, it should be converges,
but may be the unstablility of calculations.
|#

;; To check a solution is correct, just calculate: (E-G) x - b =? 0
(define (check-sol x G b)
  (vector-map - x (matrix-*-vector G x) b))
(let* ([n 3]
       ;;[G (vector #(1.2 0 0) #(0 4.3 0) #(0 0 7.4))]
       ;;[G (vector #(-2.2 0 0) #(0 -1.3 0) #(0 0 -5.4))]
       ;;对称阵,即使特征值大于1, 加速器依然可以迫使其收敛
       ;;[G  (Jordan -2.1 n)]; +/- 1.1
       ;;但是约旦标准型c>1的特征值, 加速器的压制不行, c<-1 则又可以,  why?
       [G (vector #(1.2 -1 0) #(0 -1.6 -1) #(0 0 -1.6))]
       [b (make-vector n 25.0)]
                                        ;[x0 (make-vector n 0.0)]
       [x0 (vector 1.0 2.0 3.0)])
  (newline)
  (print "----------- orginal: ")
  (for-each (lambda (i)
              (print (vector-norm (check-sol (stream-ref (multi-general-stream-from G b x0 n) i)
                                             G b))))
            (range 5))
  (print "----------- aitken: ")
  (for-each (lambda (i)
              (print (vector-norm (check-sol (stream-ref (multi-Aitken-general-strm G b x0 n) i)
                                             G b))))
            (range 5))
  (print "----------- super:  ")
  (for-each (lambda (i)
              (print (vector-norm (check-sol (stream-ref (super-multi-general-strm G b x0 n) i)
                                             G b))))
            (range 30))
  )


